/*
**  File: AquaFuncts_dr.c
**
**  This file is part of the AQUA program suite, which is being
**  developed in the NMR Spectroscopy Department, Bijvoet Center for 
**  Biomolecular Research, Utrecht. 
**
**
**  AUTHOR(S) : Ton Rullmann, Jurgen Doreleijers
**  VERSION   : 1.5.0
**  DATE      : Dec 17, 1997
**
**  This file contains routines for handling Aqua distance restraint files
**  and distance restraint data structures.
**
*/

#include <stdlib.h>
#include <stdio.h>
#include <string.h>
#include "Compiler.h"
#include "General.h"

#include "AquaMacros.h"
#include "AquaTypes.h"

/****************  function prototypes  ****************/

#include "AquaRange.h"
#include "AquaFuncts.h"
#include "AquaFuncts_io.h"
#include "AquaStrucset.h"
#include "AquaFiles.h"
#include "MenuFuncts.h"
#include "Functs.h"

#include "AquaFuncts_dr.h"

static void Swap_Atoms( DRestraint *d );
static void Sort_List( DRestraint *DRs, int DR_count );
static Boolean Restraint_In_Order( DRestraint a, DRestraint b );
static Boolean Atoms_In_Order( Atom at1, Atom at2, Boolean *same );   
static Boolean CheckDRMultipl( DRestraint *DRs, DRestraint *DRt, 
                               const int count, const int kount,
                               float *newlo, float *newup, Bndstat *newstat,
                               int *ndoubl, int *nequal );
static void CleanupDR_err1( const int count, DRestraint *DRs );
static void CleanupDR_lst1( const int count, DRestraint *DRs );
static void CleanupDR_err2( const int count, DRestraint *DRs,
                            const int kount, DRestraint *DRt );
static void CleanupDR_lst2( const int count, DRestraint *DRs,
                            const int kount, DRestraint *DRt );
static void WriteDR( FILE *fp, DRestraint *DRs );
static Bndstat CheckBoundsValTyp( DRestraint *DRt );

/****************  globals  ****************/

extern Boolean  iaf;                           /* import from Functs */
extern Card     msg_string;                    /* import from Functs */
extern int      terseness;                     /* import from Functs */
extern Range    *range;                        /* import from AquaRange */
extern char     *project_name;                 /* import from AquaFuncts */

/* BEGIN_EXTH */

/****************  InitDRList  ****************/
/*
*/
 
void InitDRList( DRList *list )
{
    list->dist_count = 0;
    list->struc_count = 0;
    list->list = NULL;
}

/****************  ProcessDRInput  ****************/
/*
**  Process distance restraint input.
*/

Boolean ProcessDRInput( DRList *DR_list, char *DR_type, Card DR_ident )
{
    char        *fct_name = "ProcessDRInput";
    int         item, count;
    Filtyp      fil_type;
    DRestraint  *DRs;
    Filnam      filnam;
    Boolean     done = false;
    FILE        *fp = NULL;
    Menu        *oldmenu;
    
    static MenuItem menulist[] =
    {
        { item_aqua,   "aqua",   none, on }, 
        { item_dgo,    "dg_old", none, off },
        { item_biosym, "biosym", none, on }, 
        { item_END,    "end",    none, on } 
    };
    static Menu menu =
    {
        "ProcessDRInput menu", "format", "Give file format:", 
        sizeof(menulist)/sizeof(MenuItem)
    };
    
    if ( terseness < 6 ) {
        STATUS( "Beginning restraint input" );       
    }
    
/* save menu and setup new one */
    
    oldmenu = SetMenu( &menu, menulist );

/* determine file type */

    while ( true )
    {
        item = ReadMenuItem();
        switch ( item )
        {
            case item_aqua:
                fil_type = aqua_drfil_type;
                goto data;
            case item_dgo:
                fil_type = dgo_drfil_type;
                goto data;
            case item_biosym:
		if ( NOTEQ( DR_type, "noe" ) )
		{
		    puts( "biosym option implemented only for noe restraints" );
		    break;
		}
                fil_type = biosym_drfil_type;
                goto data;
            case item_END:
            case item_NULL:
                goto end;                  /* no processing done */
            default:
                ProcessMenuError( menu_fatal );
                break;
        }    
    }

/* read the file name and open the file */

data:
    SetOneMenu( "Give name of distance-restraint file (or - to end):" );
    while ( !done )
    {
	MReadWord( filnam );
        if ( *filnam == '-' || strlen( filnam ) == 0 )
	    goto end;                      /* no processing done */
	if ( ! ( done = OPENr( fp, filnam ) ) )
	{
	    printf( " * Restraint-file <%s> could not be opened\n", filnam );
	    MError();
	}
    }
    if ( terseness < 6 ) {
        printf( "Restraint-file <%s> opened for reading in <%s> format\n",
            filnam, FiltypNam( fil_type ) );
    }

/* determine the number of distance restraints */

    count = DRestraintCount( fp, fil_type );
    if ( count <= 0 )
    {
        printf( " * Number of restraints determined as: %i\n", count );
	puts( " * No restraint processing" );
        if ( !iaf )
            exit( 10 );                    /* special exit status */
        goto end;                          /* no processing done */
    }
    else {
        if ( terseness < 6 ) {
            printf( "Will read %i restraints\n", count );
        }
    }

/* now read the file */

    DRs = ReadDRestraints( fp, fil_type, DR_type, count, DR_ident );
    if ( !DRs )
    {
        puts( " * Error while reading distance restraints" );
        if ( !iaf )
            exit( 1 );
        goto end;                          /* no processing done */
    }
    else {
        if ( terseness < 6 ) {
	    printf( " * Read %i distance restraints successfully!\n", count );
        }
    }
/* done */

    ResetMenu( oldmenu );                  /* reinstall the saved menu */
    if ( terseness < 6 ) {
        STATUS( "Restraint input finished" );       
    }
    DR_list->list = DRs;
    DR_list->dist_count = count;
    DR_list->struc_count = 1;
    return( true );                        /* return the list pointer and count */

end:
    ResetMenu( oldmenu );
    return( false );
}

/****************  MapDRInput  ****************/
/*
**  Map the atom names after ProcessDRInput.
*/

Boolean MapDRInput( DRList *DR_list, Mol *molptr, const int struc_count )
{
    char       *fct_name = "MapDRInput";
    char        *lib_nam_temp;
/*    DRestraint *DRs; */
/*    int        i; */

    if ( DR_list->struc_count != 1 )
	FATAL( "Can't handle multiple-model DRestraint list" );

    if ( terseness < 6 ) {
        STATUS( "mapping atom names" );
    }
    
    if ( OptionTrue( opt_atmlib ) ) {
        lib_nam_temp = "user";
    } else {
        lib_nam_temp = "";
    }
    if ( !MapAtomsInDRestraints( DR_list, ext_to_int, 
			   lib_nam_temp, 
			   molptr, struc_count ) )
    {
	ERROR( "map failed - destroying the restraint list" );
	free( DR_list->list );               /* destroy the list */
	InitDRList( DR_list ); 
	return( false );
    }

    if ( OptionTrue( opt_prtrst ) )
	PrintDRList( DR_list );

    return( true );
}

/****************  ProcessDROutput  ****************/
/*
**  Process distance restraint output.
*/

void ProcessDROutput( DRList *DR_list, const int skip_count, char *DR_type,
                      const Boolean toscreen, Mol *molptr, const int struc_count )
{
    char     *fct_name = "ProcessDROutput";
    Boolean  done = false;
    Filtyp   fil_type = no_fil_type;
    Filnam   data_filnam;
    FILE     *data_fp = NULL;
    Menu     *oldmenu;
    int      item = 0;
    char     *rstout, *lib_nam;
    Card     ident;

    static MenuItem menulist[] =
    {
        { item_aqua,    "aqua",    none, on }, 
        { item_disman,  "disman",  none, off }, 
        { item_disgeo,  "disgeo",  none, off }, 
        { item_diana,   "diana",   none, on }, 
        { item_xplor,   "xplor",   none, on }, 
        { item_biosym,  "biosym",  none, on }, 
        { item_nmrstar, "nmrstar", none, on }, 
        { item_END,     "end",     none, on }  
    };                      /* same as in AquaWhat, AquaHow, AquaPseudo and GetLIBName */
    static Menu menu =
    {
        "ProcessDROutput menu", "format", "Give file format:", 
        sizeof(menulist)/sizeof(MenuItem)
    };
    

    if ( terseness < 6 ) {
        STATUS( "Beginning restraint output" );       
    }

/* save menu and setup new one */
    
    oldmenu = SetMenu( &menu, menulist );
    
/* select the output format */

    while ( !done )
    {
        item = ReadMenuItem();
        switch ( item )
        {
            case item_aqua:
                fil_type = aqua_drfil_type;
                done = true;
		break;
            case item_diana:
                fil_type = diana_drfil_type;
                done = true;
		break;
            case item_xplor:
                fil_type = xplor_drfil_type;
                done = true;
		break;
            case item_biosym:
                fil_type = biosym_drfil_type;
                done = true;
		break;
            case item_nmrstar:
                fil_type = nmrstar_drfil_type;
                done = true;
		break;
            case item_END:
            case item_NULL:
                ResetMenu( oldmenu );
                return;                  /* no processing done */
            default:
                ProcessMenuError( menu_fatal );
                break;
        }    
    }
    rstout = GetMenuItemKey( item );

    if ( item == item_aqua || item == item_nmrstar )
    {
	if ( iaf )
	    puts( "Now choose the atom name format:" );
	rstout = GetLIBName( item_NONE );
	if ( !rstout )
	    goto end;
    }

/* read the file name and open the data file */

    SetOneMenu( "Give name of data file (or - to end):" );
    done = false;
    if ( toscreen )
    {
	data_fp = stdout;
	done = true;
    }

    while ( !done )
    {
	MReadWord( data_filnam );
        if ( *data_filnam == '-' || strlen( data_filnam ) == 0 )
	    goto end;                       /* no processing */

	if ( ( data_fp = fopen( data_filnam, "w" ) ) == NULL )
	{
	    printf( " * Data-file <%s> could not be opened\n", data_filnam );
	    MError();
	}
	else
	    done = true;
    }

/* first map, then write the data */

    lib_nam = rstout;
    if ( EQ( rstout, "aqua" ) && fil_type == aqua_drfil_type )
	lib_nam = "copy";
    if ( !MapAtomsInDRestraints( DR_list, int_to_ext, lib_nam, 
                                 molptr, struc_count ) )
    {
	printf( " * Skipping restraint output\n" );
	goto end;
    }

    SetOneMenu( "Give text for title line (between \''s):" );
    MReadText2( ident, size_LINE, '\'' );
    if ( WriteDRestraints( data_fp, fil_type, DR_list, DR_type, skip_count, 
                           ident, rstout ) ) {
        if ( terseness < 6 ) {
	    printf( " * Distance restraints have been written\n" );
        }
    } else {
	printf( " * Distance restraints have not been written\n" );
    }
    
end:
    if ( data_fp && data_fp != stdout )
	fclose( data_fp );
    ResetMenu( oldmenu );                    /* reinstall the saved menu */
    if ( terseness < 6 ) {
        STATUS( "Restraint output finished" );       
    }
}

/****************  DRestraintCount  ****************/
/*
**  Determine the number of distance restraints from a data file.
*/

int DRestraintCount( FILE *fp, const Filtyp type )
{
    int    count;
    
    switch ( type )
    {
	case aqua_drfil_type:
	    FlushLine( fp );
	    fscanf( fp, "count %i", &count );
	    rewind( fp );
	    break;
        case biosym_drfil_type:
            count = NoeCount_biosym( fp );
            break;
        default:
	    printf( " * Reading of distance restraint file type <%s> not implemented\n", 
	            FiltypNam( type ) );
	    exit( 1 );
    }
    
    return( count );
}

/****************  ReadDRestraints  ****************/
/*
**  Read distance restraints from a data file.
*/

DRestraint *ReadDRestraints( FILE *fp, const Filtyp type, char *DR_type,
                             const int count, Card ident )
{
    int         DR_read;
    DRestraint  *DRs;
    
    DRs = CREATE_NEW( DRestraint, count );
    InitDRestraints( DRs, count );
    
    switch ( type )
    {
        case aqua_drfil_type:
            DR_read = ReadRestraintFile( fp, DRs, count, ident, DR_type, GetOptionIVal( opt_shift ) );
            break;
        case biosym_drfil_type:
            DR_read = ReadNoes_biosym( fp, DRs, count, ident );
            break;
	default:
	    printf( " * Reading of distance restraint file type <%s> not implemented\n", 
	            FiltypNam( type ) );
	    exit( 1 );
    }
    fclose( fp );
    
    if ( DR_read != count )
    {
        printf( " * Number of restraints read (%i) does not match the expected number (%i)",
                DR_read, count );
	free( DRs );               /* destroy the list */
        return( NULL );
    }

    return( DRs );
}

/****************  WriteDRestraints  ****************/
/*
**  Write distance restraints to a data file.
*/

Boolean WriteDRestraints( FILE *fp, const Filtyp type, DRList *theList,
                          char *DR_type, const int skip_count,
		          Card ident, char *format )
{
    char       *fct_name = "WriteDRestraints";
    char       *Mol_label = GetVar( "AQUAPROJECT", false );
    Boolean    ok = true;
    DRestraint *DR_all;
    int        DR_count;

    if ( theList->struc_count != 1 )
	FATAL( "Can't handle multiple-model DRestraint list" );
	
    DR_all = theList->list;
    DR_count = theList->dist_count;
    

    switch( type )
    {
	case aqua_drfil_type:
	    WriteRestraintFile( fp, DR_all, DR_count, skip_count, ident, DR_type, format );
	    break;
	case diana_drfil_type:
	    WriteRestraints_diana( fp, DR_all, DR_count, skip_count, ident );
	    break;
	case xplor_drfil_type:
	    WriteRestraints_xplor( fp, DR_all, DR_count, skip_count, ident );
	    break;
	case biosym_drfil_type:
	    WriteRestraints_biosym( fp, DR_all, DR_count, skip_count, ident, DR_type );
	    break;
	case nmrstar_drfil_type:
	    ok = WriteRestraints_nmrstar( fp, DR_all, DR_count, skip_count,  
	                                  ident, DR_type, Mol_label );
	    break;
	default:
	    printf( " * Writing of distance restraint file type <%s> not implemented\n", 
	            FiltypNam( type ) );
	    exit( 1 );
    }
    return( ok );
}

/****************  InitDRestraints  ****************/
/*
**  Initialize distance restraints.
*/

void InitDRestraints( DRestraint *DRs, int count )
{
    while ( count-- )
    {
	InitAtom( &(DRs->atom_A) );
	InitAtom( &(DRs->atom_B) );
	DRs->lo_bnd = 0.0;
	DRs->up_bnd = 0.0;
	DRs->dist = 0.0;
	DRs->keep = true;
	DRs++;
    }
}

/****************  PrintDRList  ****************/
/*
**  Print a DRList.
*/

void PrintDRList( DRList *DR_list )
{
    DRestraint *DRs;
    int        i;
/*    char       *typ; */

    DRs = DR_list->list;
    for ( i = 1; i <= DR_list->struc_count; i++ )
    {
	printf( "Model %i\n", i );
	PrintDRestraints( DRs, DR_list->dist_count );
	DRs += DR_list->dist_count;
    }
}

/****************  PrintDRestraints  ****************/
/*
**  Print a list of DRestraints.
*/

void PrintDRestraints( DRestraint *DRs, int count )
{
    char   *typ;

    while ( count-- )
    {
	switch ( DRs->stat )
	{
	    case bnd_upper:
		typ = "upp";
		break;
	    case bnd_lower:
		typ = "low";
		break;
	    case bnd_uplo:
		typ = "uplo";
		break;
	    case bnd_uld:
		typ = "uld";
		break;
	    case bnd_ok:
		typ = "ok";
		break;
	    default:
		typ = "err";
	}
	printf( "%-4s ", typ );
	WriteFAtom( stdout, &DRs->atom_A, use_both );
	putchar( ' ' );
	WriteFAtom( stdout, &DRs->atom_B, use_both );
	printf( " %f %f %f\n", DRs->lo_bnd, DRs->up_bnd, DRs->dist );
        DRs++;
    }
}

/****************  MapCopyInDRestraints  ****************/
/*
**  Copy external to internal names.
*/

void MapCopyInDRestraints( DRestraint *DRs, int count )
{
    while ( count-- )
    {
        memcpy( DRs->atom_A.atom_int, DRs->atom_A.atom_ext, size_ID );
        memcpy( DRs->atom_B.atom_int, DRs->atom_B.atom_ext, size_ID );
        DRs++;
    }
}

/****************  MapAtomsInDRestraints  ****************/
/*
**  For a DRestraint list map the 'external' atom names to the 
**  'internal' (Aqua) ones, or vice versa (determined by <how>).
**  If needed, multiple lists (one per model) are created and a loop 
**  over them is performed.
**  Cf. MapAtomsInSet.
*/

Boolean MapAtomsInDRestraints( DRList *theList, const int how, char *lib_nam, 
                               Mol *molptr, const int struc_count )
{
    char        *fct_name = "MapAtomsInDRestraints";
    Atomptr     *atomptrs;
    DRestraint  *newDR, *newDRs, *DRs;
    Filnam      predef_nam;
    int         DRcount, count;

    if ( theList->struc_count != 1 )
	FATAL( "Can't handle multiple-model DRestraint list" );

    DRs = theList->list;
    DRcount = theList->dist_count;

    if ( terseness < 6 ) {
        STATUS( "preparing to map atom names in the restraint list" );
    }
    
/* allocate a list of pointers to the atoms of the restraint list */

    atomptrs = CREATE_NEW( Atomptr, 2*DRcount );
    if ( !atomptrs )
	NOMEM( "atomptrs" );

/* predefine the library name for use with the "whatif" option */

    strcpy( predef_nam, project_name );
    strcat( predef_nam, FiltypExt( atlib_type ) );

/* DO THIS HERE ? >>>>>>>>>>>>>>>>>>>>>>>> */

    if ( EQn( lib_nam, "whatif", 6 ) && struc_count > 1 )

/* use the model-dependent molecules */
    {
	newDRs = CREATE_NEW( DRestraint, DRcount * struc_count );
	if ( !newDRs )
	    NOMEM( "newDRs" );
	newDR = newDRs;
	count = 0;
	while ( ++count <= struc_count )
	{
	    memcpy( newDR, DRs, (size_t)DRcount * sizeof( DRestraint ) );
	    sprintf( predef_nam, "%s_%03i%s", project_name, 
		     count, FiltypExt( atlib_type ) );
	    if ( !MapAtomsIn1DR( newDR, DRcount, atomptrs, how, lib_nam, 
	                         predef_nam, molptr ) )
		return( false );
	    newDR += DRcount;
	}

	free( DRs );                   /* save the new list */
	theList->list = newDRs;
	theList->dist_count = DRcount;
	theList->struc_count = struc_count;
    }

/* or use the unique molecule */

    else
    {
	if ( !MapAtomsIn1DR( DRs, DRcount, atomptrs, how, lib_nam, predef_nam, molptr ) )
	    return( false );
    }

    return( true );
}

/****************  MapAtomsIn1DR  ****************/
/*
**  For a DRestraint list map the 'external' atom names to the 
**  'internal' (Aqua) ones, or vice versa (determined by <how>).
**  Cf. MapAtomsInMol.
**
**  The library name <lib_nam> is derived from the restraint file type <rstinp>, 
**  if it is "". Names are just copied if <lib_nam> is "copy" or "aqua". 
**
**  <molptr> is needed to expand residue wildcards in the AtomLIB and for
**  setting the residue name if it is "UNK", but can be NULL if <lib_nam> is
**  "copy" or "aqua".
**
**  Name guessing - as implemented in MapPtrAtoms->FindAtom - is tried 
**  if <opt_guess> has been set.
*/

Boolean MapAtomsIn1DR( DRestraint *DRall, const int DRcount, Atomptr *atomptrs, 
                       const int how, char *lib_nam_orig, Filnam predef_nam, 
		       Mol *molptr )
{
    char        *fct_name = "MapAtomsIn1DR";
    Atomptr     *list;
    DRestraint  *DRs;
    char        *res_nam, *res_id;
    int         count;
    Mol         *atmap_mol = NULL, *atmap_mol2 = NULL;
    FctFindAt   Find = NULL;
    Boolean     guess;
    Strucset    tmpSet;
    Card        lib_nam;
    char        *lib_nam1 = lib_nam, *lib_nam2;
    
/* make a list of pointers to the atoms of the restraint list */

    list = atomptrs;
    DRs = DRall;
    count = DRcount;
    while ( count-- )
    {
	*list++ = &(DRs->atom_A);
	*list++ = &(DRs->atom_B);
        DRs++;
    }

/* make a list of atom name pairs from a library file */
/* the library name is derived from the file type, if it is "" */ 
/* if the file type is "aqua", names will just be copied */
/* else, first correct the residue names if they are "UNK" */

/* was a combination of libraries specified? if so, split the name on the + character */

    strcpy( lib_nam, lib_nam_orig );
    if ( ( lib_nam2 = strchr( lib_nam, '+' ) ) )
    {
	*lib_nam2 = '\0';      /* this sets lib_nam1 to a correct string */
	lib_nam2++;
    }

    if ( NOTEQ( lib_nam1, "copy" ) && NOTEQ( lib_nam1, "aqua" ) )
    {
	if ( EQ( lib_nam1, "" ) )
	{
	    if ( !( lib_nam1 = GetLIBName( opt_rstinp ) ) )
		return( false );
	}
	if ( NOTEQ( lib_nam1, "copy" ) && NOTEQ( lib_nam1, "aqua" ) )
	{
	    if ( !molptr )
	    {
		ERROR( "no molecule defined" );
		return( false );
	    }
	    atmap_mol = ProcessAtomLIB( molptr, lib_nam1, predef_nam );
    
	    DRs = DRall;
	    count = DRcount;
	    while ( count-- )
	    {
		res_nam = DRs->atom_A.res_nam;
		res_id  = DRs->atom_A.res_id;
		if ( EQID( res_nam, "UNK  " ) )
		{
		    memcpy( res_nam, FindResidueNamFromId( res_id, molptr ), size_ID );
		    MakeUniqResid( &DRs->atom_A );
		}
		res_nam = DRs->atom_B.res_nam;
		res_id  = DRs->atom_B.res_id;
		if ( EQID( res_nam, "UNK  " ) )
		{
		    memcpy( res_nam, FindResidueNamFromId( res_id, molptr ), size_ID );
		    MakeUniqResid( &DRs->atom_B );
		}
		DRs++;
	    }
	}
    }

/* second library - handling same as in MapAtomsInMol */

    if ( lib_nam2 && atmap_mol )
    {
	if ( IsEmptyStr( lib_nam2 ) )
/*	    lib_nam2 = FiltypLIBNam( molptr->inp_filtyp );   */
	    lib_nam2 = GetLIBName( item_NONE );              
	if ( ( atmap_mol2 = ProcessAtomLIB( molptr, lib_nam2, predef_nam ) ) )
	{
	    StrucsetInit( &tmpSet, "allDefs", false );
	    StrucsetCreate( &tmpSet, atmap_mol->atom_count + atmap_mol2->atom_count, 0 );
	    AddAtomSToSet( &tmpSet, atmap_mol->atoms, atmap_mol->atom_count );
	    AddAtomSToSet( &tmpSet, atmap_mol2->atoms, atmap_mol2->atom_count );
	    DeleteMol( atmap_mol2 );
	    if ( OptionTrue( opt_find2 ) )
		MakeResidues( &tmpSet, false );     /* cf. MakeDefsMol */
	    atmap_mol = tmpSet.molecule;
	    sprintf( msg_string, "%i atom name definitions active for present molecule", 
		 atmap_mol->atom_count );
            if ( terseness < 6 ) {
	        STATUSm;
            }
	    if ( OptionTrue( opt_prtlib ) )
		WriteFAtoms( stdout, atmap_mol->atoms, atmap_mol->atom_count );
	}
    }

/* print what is going to happen */

    if ( !atmap_mol )
    {
        if ( terseness < 6 ) {
	    if ( how == int_to_ext )
	        STATUS( "copying internal to external atom names" );
	    else
	        STATUS( "copying external to internal atom names" );
        }
    }
    else
    {
	if ( how == int_to_ext )
	{
            if ( terseness < 6 ) {
	        STATUS( "mapping internal to external atom names" );
            }
	    if ( OptionTrue( opt_find2 ) )
	    {
                if ( terseness < 6 ) {
		    STATUS( "using routine FindAtom2" );
                }
		Find = FindAtom2;
	    }
	    else
	    {
                if ( terseness < 6 ) {
		    STATUS( "using routine FindAtom3" );
                }
		Find = FindAtom3;
	    }
	}
	else
	{
            if ( terseness < 6 ) {
	        STATUS( "mapping external to internal atom names" );
            }
	    if ( OptionTrue( opt_find2 ) )
	    {
                if ( terseness < 6 ) {
		    STATUS( "using routine FindAtom2" );
                }
		Find = FindAtom2;
	    }
	    else
	    {
                if ( terseness < 6 ) {
		    STATUS( "using routine FindAtom3" );
                }
		Find = FindAtom3;
	    }
	}
    }
    guess = OptionTrue( opt_guess ) && how == ext_to_int;
    if ( guess )
    {
        if ( terseness < 6 ) {
	    STATUS( "guessing ON" );
        }
	ModOption( opt_rstinp, lib_nam1 );       /* save for use in GuessAtom */
    }
    else {
        if ( terseness < 6 ) {
	    STATUS( "guessing OFF" );
        }
    }
    
/* now do the mapping or copying */

    MapPtrAtoms( atomptrs, 2*DRcount, atmap_mol, how, guess, Find );

/* free the list again */

    if ( atmap_mol )
	DeleteMol( atmap_mol );

    return( true );
}

/****************  MakeDRList  ****************/
/*
**  Create a DRList from a Distinf/Distrst list-pair.
*/

void MakeDRList( DRList *theList, Distinf *dist_inf, Distrst *dist_rst, 
                 const int dist_count, Atom *atoms )
{
    char        *fct_name = "MakeDRList";
    DRestraint  *DR_all, *DRs;
    int         count;

    DR_all = CREATE_NEW( DRestraint, dist_count );
    if ( !DR_all )
	NOMEM( "DR_all" );
    
    DRs = DR_all;
    count = dist_count;
    while ( count-- )
    {
	memcpy( &(DRs->atom_A), &atoms[dist_inf->idx_A], sizeof( Atom) );
	memcpy( &(DRs->atom_B), &atoms[dist_inf->idx_B], sizeof( Atom) );
	DRs->lo_bnd = dist_rst->lo_bnd;
	DRs->up_bnd = dist_rst->up_bnd;
	DRs->dist   = dist_rst->dist;
	DRs->idx_A = dist_inf->idx_A;
	DRs->idx_B = dist_inf->idx_B;
	DRs->keep = true;
	switch( dist_rst->typ_bnd )
	{
	    case 'u':
		DRs->stat = bnd_upper;
		break;
	    case 'l':
		DRs->stat = bnd_lower;
		break;
	    case 'b':
		DRs->stat = bnd_uplo;
		break;
	    case 'B':
		DRs->stat = bnd_uld;
		break;
	    default:
		printf( "* unexpected bounds type <%c>\n", dist_rst->typ_bnd );
		exit( 1 );
	}
	dist_inf++;
	dist_rst++;
        DRs++;
    }

    theList->list = DR_all;
    theList->dist_count = dist_count;
    theList->struc_count = 1;
}

/****************  SelectRestraints  ****************/
/*  
**  Select restraints according to a Range definition.
**  The total number of skipped restraints in the final list
**  (i.e. those skipped before, plus new skips) is returned.
**  The Range is a global variable, used in routine InResidueRange.
*/

int SelectRestraints( DRList *DR_list, Mol *molptr )
{
    char       *fct_name = "SelectRestraints";
    DRestraint *DRs = DR_list->list;
    int        DR_count = DR_list->dist_count;
    Menu       *oldmenu;
    FILE       *fp = NULL;
    Filnam     filnam;
    Boolean    done = false;
    int        count, new_count = 0, idx_A, idx_B, how_sel;
    Boolean    select_A, select_B;
    Residue    *residues = molptr->residues;

    if ( DR_list->struc_count != 1 )
	FATAL( "Can't handle multiple-model DRestraint list" );

    range = NULL;
    oldmenu = SetOneMenu( "Give name of ranges file (or - to end):" );
    while ( !done )
    {
	MReadWord( filnam );
	if ( *filnam == '-' || strlen( filnam ) == 0 )
	    break;
	if ( ! ( done = OPENr( fp, filnam ) ) )
	{
	    printf( " * ranges file <%s> could not be opened\n", filnam );
	    MError();
	}
    }
    ResetMenu( oldmenu );
    if ( done )
    {
	fclose( fp );
	range = ReadRanges( filnam );
    }
    if ( !range )
    {
	count = DR_count;
	while ( count-- )
	{
	    if ( DRs->keep )
		new_count++;
	    DRs++;
	}
	return( DR_count - new_count );
    }

/* how to select restraints? */

    how_sel = GetOptionSubItem( opt_rest_sel );
    switch ( how_sel )
    {
	case restsel_both:
            if ( terseness < 6 ) {
	        puts( " * Selecting restraints for which BOTH residues are in allowed range" );
            }
	    break;
	case restsel_one:
            if ( terseness < 6 ) {
	        puts( " * Selecting restraints for which AT LEAST ONE residues is in allowed range" );
            }
	    break;
	default:
	    FATAL( "Illegal opt_rest_sel value - program error" );
    }

/* cf. mechanism in RangeDistLists */

    count = DR_count;
    while ( count-- )
    {
	if ( !DRs->keep )
	{
	    DRs++;
	    continue;
	}
/* no mapping done yet; so use the original residue number and chain
   to find the residue */
	idx_A = FindResidueIndexFromOrigNum( DRs->atom_A.res_num_orig, 
	                  DRs->atom_A.chain, DRs->atom_A.res_insert, molptr );
	idx_B = FindResidueIndexFromOrigNum( DRs->atom_B.res_num_orig, 
	                  DRs->atom_B.chain, DRs->atom_B.res_insert, molptr );
	if ( idx_A < 0 || idx_B < 0 )
	{
	    printf( "Residue(s) not found: " );
	    PrintDRestraints( DRs, 1 );
	}
	select_A = idx_A < 0 ? false : InResidueRange( residues+idx_A );
	select_B = idx_B < 0 ? false : InResidueRange( residues+idx_B );
	if (  how_sel == restsel_both ? ( select_A && select_B ) : ( select_A || select_B ) )
	{
	    new_count++;
	    DRs->keep = true;
	}
	else
	    DRs->keep = false;
	DRs++;
    }
    
    if ( terseness < 6 ) {
        printf( " * Keeping %i distance restraints\n", new_count );
    }
    return( DR_count - new_count );
}

/****************  SortRestraints  ****************/
/*  
**  Sorts restraint list in ascending order.
**  If <classify> the list is also partitioned in FIVE separate lists:
**  intra (i=j), sequential or short (|i-j|=1), 
**  medium (1<|i-j|<=max_MEDIUMRANGE), long (|i-j|>max_MEDIUMRANGE) 
**  and inter-chain.
**
**  AUTHOR: J.F. Doreleijers. Modified by Ton Rullmann.
*/

void SortRestraints( DRList *DR_list, Boolean classify )
{
    char           *fct_name = "SortRestraints";
    DRestraint     *DRs = DR_list->list;
    int            DR_count = DR_list->dist_count;
    int            num_range;
    int	           i, j, k, sum, diff_resnum;
    int	           *list_count;
    Boolean	   same, insert_present = false;
    ContactRange   contact_range;
    DRestraint     **DRl;


    if ( DR_list->struc_count != 1 )
	FATAL( "Can't handle multiple-model DRestraint list" );

    num_range = classify ? num_CONTACTRANGE : 1;

    DRl        = CREATE_NEW( DRestraint *, num_range );
    list_count = CREATE_NEW( int, num_range );
    if ( !list_count || !DRl )
	NOMEM( "lists" );

    /* eg DRl[1][123] is the 124 st restraint in the second sequential list */
     
    for ( i=0; i<num_range; i++ ) {
	list_count[i] = 0;
	DRl[i] = CREATE_NEW( DRestraint, DR_count );   /* Allocate second dimension */
	if ( !DRl[i] )
	    NOMEM( "DRl" );
    }
			      
    if ( classify )
    {
        if ( terseness < 9 ) {
	    STATUS( "Start partitioning restraint list" );
        }
    
	for ( i=0; i<DR_count; i++ ) {
	    
	    if ( DRs[i].atom_A.res_insert != ' ' || DRs[i].atom_B.res_insert != ' ' )
		insert_present = true;
    
	    diff_resnum = DRs[i].atom_A.res_num_orig - DRs[i].atom_B.res_num_orig ;
	    
	    if ( diff_resnum < 0 )
		diff_resnum *= -1; 
		    
	    contact_range = undef_rng;
	    if ( diff_resnum == 0 )
		contact_range = intra_rng;
	    if ( diff_resnum == 1 )
		contact_range = seq_rng;
	    if ( ( diff_resnum > 1 ) && ( diff_resnum <= max_MEDIUMRANGE ) )
		contact_range = medium_rng;
	    if ( diff_resnum > max_MEDIUMRANGE )
		contact_range = long_rng;
	    if ( NOTEQID( DRs[i].atom_A.chain, DRs[i].atom_B.chain ) )
		contact_range = inter_rng;
	    
	    switch ( contact_range ) {
	    
		case intra_rng :
		    DRl[0][ list_count[0] ] = DRs[i];
		    ++list_count[0];
		    break;
		case seq_rng :
		    DRl[1][ list_count[1] ] = DRs[i];
		    ++list_count[1];
		    break;
		case medium_rng :
		    DRl[2][ list_count[2] ] = DRs[i];
		    ++list_count[2];
		    break;
		case long_rng :
		    DRl[3][ list_count[3] ] = DRs[i];
		    ++list_count[3];
		    break;
		case inter_rng :
		    DRl[4][ list_count[4] ] = DRs[i];
		    ++list_count[4];
		    break;
		default :
		    FATAL( "undefined contact_range" );
	    }
	}
    }
    else
    {
	for ( i=0; i<DR_count; i++ ) {
	    DRl[0][ list_count[0] ] = DRs[i];
	    ++list_count[0];
	}
    }
    
    sum = 0;
    
    for ( i=0; i<num_range; i++) 
	sum += list_count[i];
	
    if ( sum != DR_count )
	FATAL( "bad sum" );
	
    k = 0;
    
    for ( i=0; i<num_range; i++ ) {                    /* Loop over the lists */
    
	sprintf( msg_string, "Sorting list %d with %d restraints", i, list_count[i] );
        if ( terseness < 6 ) {
	    STATUSm;
        }
	
	for ( j=0; j<list_count[i]; j++ )      /* Swap atom A and B if necessary */
	{
	    if ( ! Atoms_In_Order( DRl[i][j].atom_A, DRl[i][j].atom_B, &same ) )
		    Swap_Atoms( &(DRl[i][j]) );
	}
	Sort_List( DRl[i], list_count[i] );             /* Sort the lists */
	
	for ( j=0; j<list_count[i]; j++ ) {             /* And copy them back */
	    DRs[k] = DRl[i][j];
	    k++;
	}
	
    }
    
    if ( k != DR_count )
	FATAL( "bad number" );

    for ( i=0; i<num_range; i++ )
	free( DRl[i] );
    free( DRl );
    free( list_count );
    
    if ( insert_present )
	ERROR( "sorting partially incorrect due to presence of residue insertion code(s)" );
}

/****************  CleanupDRList  ****************/
/*
**  Cleanup of DRList: 
**  - check whether atoms in restraint are different
**  - check for multiple entries of same distance
**  - combine upper and lower bounds belonging to same distance
**  The total number of skipped restraints in the final list
**  (i.e. those skipped before, plus new skips) is returned.
**  This count is over the first restraint model only
**  (cf. theList->dist_count, which also refers to one model). 
**
**  All restraint models should result in the same (apart from name swaps)
**  cleaned up list.
**
**  Special exit flag value (cf. AquaHow):
**    20: not all restraints were handled because some atoms were not found 
*/

int CleanupDRList( DRList *theList, const int skip_count, int *exit_flag )
{
    char        *fct_name = "CleanupDRList";
    DRestraint  *DRall = theList->list;
    int         struc_num, DRcount = theList->dist_count;
    int         count, kount;
    int         nerr     = 0, nnoat     = 0, nsame     = 0;
    int         ndoubl   = 0, nequal    = 0, nincomp   = 0;
    int         nerr_old = 0, nnoat_old = 0, nsame_old = 0;
    Bndstat     stat, newstat;
    float       newup, newlo;
    DRestraint  *DRs, *DRt;

    if ( terseness < 6 ) {
        STATUS( "Checking restraints" );
    }

/* initialize to keep all restraints except if an atom index is -1 */
/* and check the values */
/* we don't look at the keep field */

    DRs = DRall;
    for ( struc_num = 1; struc_num <= theList->struc_count; struc_num++ )
    {    
	if ( theList->struc_count > 1 )
	    printf( " Restraint model %i\n", struc_num );
	nerr = 0;
	nnoat = 0;
	nsame = 0;
	count = 0;
	while ( ++count <= DRcount )
	{
	    stat = CheckBoundsValTyp( DRs );
            
	    switch ( stat )
	    {
		case bnd_sml_lo:
		    printf( " * Error: lower bound smaller than %f\n", small_BOUND );
		    CleanupDR_err1( count, DRs );
		    nerr++;
		    DRs->keep = false;
		    break;
		case bnd_big_up:
		    printf( " * Error: upper bound greater than %f\n", large_BOUND );
		    CleanupDR_err1( count, DRs );
		    nerr++;
		    DRs->keep = false;
		    break;
		case bnd_incomp:
		    printf( " * Error: upper bound < lower bound\n" );
		    CleanupDR_err1( count, DRs );
		    nerr++;
		    DRs->keep = false;
		    break;
		case bnd_incompd:
		    printf( " * Error: distance outside bounds\n" );
		    CleanupDR_err1( count, DRs );
		    nerr++;
		    DRs->keep = false;
		    break;
		case bnd_noatom:
		    printf( " Warning: bound skipped because atom(s) could not be found\n" );
		    CleanupDR_lst1( count, DRs );
		    nnoat++;
		    DRs->keep = false;
		    break;
		case bnd_sameat:
		    printf( " * Error: atoms are identical\n" );
		    CleanupDR_err1( count, DRs );
		    nsame++;
		    DRs->keep = false;
		    break;
		default:
		    break;
	    }	    
	    DRs->stat = stat;
	    DRs++;
	}
	if ( struc_num > 1 )
	{
	    if ( nerr_old != nerr )
		FATAL( "Number of errors (nerr) for this restraint model differs from previous model" );
	    if ( nnoat_old != nnoat )
		FATAL( "Number of errors (nnoat) for this restraint model differs from previous model" );
	    if ( nsame_old != nsame )
		FATAL( "Number of errors (nsame) for this restraint model differs from previous model" );
	}
	nerr_old = nerr;
	nnoat_old = nnoat;
	nsame_old = nsame;
    }

    if ( theList->struc_count > 1 && OptionTrue( opt_multi ) )
    {
	STATUS( "Multiple-model DRestraint list - will ignore opt_multi" );
	ModOption( opt_multi, "off" );
    }

    if ( OptionTrue( opt_multi ) )
    {

/* check for multiple entries - only for kept entries */

	count = 0;
	DRs = DRall;
	while ( ++count <= DRcount )
	{
	    if ( !DRs->keep )
	    {
		DRs++;
		continue;
	    }
    
	    kount = 0;
	    DRt = DRall;
	    while ( ++kount < count )
	    {
		if ( !DRt->keep )
		{
		    DRt++;
		    continue;
		}
		if (
		      (( DRt->idx_A == DRs->idx_A ) && ( DRt->idx_B == DRs->idx_B ))
		   ||
		      (( DRt->idx_A == DRs->idx_B ) && ( DRt->idx_B == DRs->idx_A ))
		   )
                                                  /* atom pairs are identical */
		{
		    if ( CheckDRMultipl( DRs, DRt, count, kount, &newlo, &newup, &newstat,
					 &ndoubl, &nequal ) )
		    {
			DRt->keep = false;        /* remove DRt. not DRs */
			DRs->stat = newstat;
			DRs->lo_bnd = newlo;
			DRs->up_bnd = newup;
                           /* the distance of each restraint remains the same 
			      but <newstat> was set to bnd_uplo */
		    }
		    else
			nincomp++;
		}
		DRt++;
	    }
	    DRs++;
	}
    }


    if ( nnoat )
    {
	printf( " * Warning: %i restraints skipped because atom(s) could not be found\n", nnoat );
	*exit_flag = 20;
    }
    if ( nerr )
	printf( " * Warning: %i restraints skipped because of error(s) in value(s)\n", nerr );
    if ( nsame )
	printf( " * Warning: %i restraints skipped because the two atoms are identical\n", nsame );
    if ( nincomp )
	printf( " * Warning: %i pairs of restraints had incompatible bounds\n", nincomp );
    if ( !nnoat && !nerr && !nsame && !nincomp ) {
        if ( terseness < 9 ) {
            puts( " * All OK!" );
        }
    }


    if ( OptionTrue( opt_multi ) )
    {
	if ( !ndoubl && !nequal && !nincomp ) {
            if ( terseness < 9 ) {
	        puts( " * All restraints are unique" );
            }
        } else {
	    if ( ndoubl )
		printf( " * %i Restraints collapsed into %i because "\
			"upper and lower bounds were\n *    specified separately\n", 
			2*ndoubl, ndoubl );
	    if ( nequal )
		printf( " * %i Restraints merged with identical restraints\n", 
			nequal );
	}
    }

/* recalculate skip_count - FOR FIRST RESTRAINT MODEL ONLY! */

    count = 0;
    kount = 0;
    DRs = DRall;
    while ( ++count <= DRcount )
    {
	if ( !DRs->keep )
	    kount++;
	DRs++;
    }
    if ( kount != ( nnoat + ndoubl + nequal + nerr + nsame + skip_count ) )
	FATAL( "Counting error" );
    if ( terseness < 6 ) {
        printf( " * Saving %i distance restraints\n", DRcount - kount );
    }

    return( kount );
}


/* END_EXTH ----------------------- only private functions follow */

/****************  Swap_Atoms  ****************/
/*
**
*/

static void Swap_Atoms( DRestraint *d )
{
    Atom  atemp = d->atom_A;
    int   itemp = d->idx_A;
    d->atom_A   = d->atom_B;
    d->atom_B   = atemp;
    d->idx_A    = d->idx_B;
    d->idx_B    = itemp;
}

/****************  Sort_List  ****************/
/*
**  Code as in "piksrt.c", described in "Numerical Recipes in C", 
**  W.H. Press et al., Cambridge University Press.
**  (Yes, it's not the most efficient sort routine possible.)
*/

static void Sort_List( DRestraint *DRs, int DR_count )
{
    int i, j;
    DRestraint temp;
    
    for ( j=1; j<=DR_count-1; j++ ) {
	temp = DRs[j];
	i = j-1;
	while ( ( i >= 0 ) && Restraint_In_Order( temp, DRs[i] ) ) {
	    DRs[i+1] = DRs[i];
	    i--;
	}
	DRs[i+1] = temp;
    }
}

/****************  Restraint_In_Order  ****************/
/*  
**  Returns true if first restraint is already before second one;
**  returns false if they should be swapped and NO LONGER 
**  exits with error if they are the same.
**  If they are the same the routine returns true.
*/

static Boolean Restraint_In_Order( DRestraint a, DRestraint b )
{
    Boolean state, same;
    
    state = Atoms_In_Order(a.atom_A, b.atom_A, &same);

    if ( !same ) 
	return ( state );

    return ( Atoms_In_Order(a.atom_B, b.atom_B, &same) );
}

/****************  Boolean Atoms_In_Order  ****************/
/* 
**  Returns true if first atom is already before second 
**  and false if they should be swapped.
*/

static Boolean Atoms_In_Order( Atom at1, Atom at2, Boolean *same )    
{         
    *same = false;
    
    if ( memcmp ( at1.chain, at2.chain, size_ID ) < 0 )
	return ( true );
    if ( memcmp ( at1.chain, at2.chain, size_ID ) > 0 )
	return ( false );
	
    if ( at1.res_num_orig < at2.res_num_orig ) 
	return ( true );
    if ( at1.res_num_orig > at2.res_num_orig ) 
	return ( false );

    if ( memcmp ( at1.atom_int, at2.atom_int, size_ID ) < 0 )
	return ( true );
    if ( memcmp ( at1.atom_int, at2.atom_int, size_ID ) > 0 )
	return ( false );
    else
	*same = true; /* chain, residue number and atom name are all the same */
   
    return ( true );
}

/****************  CheckDRMultipl  ****************/
/*
**  
*/

static Boolean CheckDRMultipl( DRestraint *DRs, DRestraint *DRt, 
                               const int count, const int kount, 
                               float *newlo, float *newup, Bndstat *newstat, 
			       int *ndoubl, int *nequal )
{
    char        *fct_name = "CheckDRMultipl";
    enum ck { ck_no, ck_same, ck_dif };
    enum ck     lochk = ck_no, upchk = ck_no;
    float       loval, upval;
    Boolean     lo, up;

/* check stat values for legality */

    switch ( DRs->stat )
    {
	case bnd_upper:
	case bnd_lower:
	case bnd_uplo:
	case bnd_uld:
	    break;
	default:
	    FATAL( "Program error - illegal value DRs->stat" );
    }
    switch ( DRt->stat )
    {
	case bnd_upper:
	case bnd_lower:
	case bnd_uplo:
	case bnd_uld:
	    break;
	default:
	    FATAL( "Program error - illegal value DRt->stat" );
    }

/* lower bound check */

    lo = true;                 /* most of the time we have a lower bound */
    loval = -99.9; /* Just for compiler warning */
    if (   ( DRs->stat == bnd_lower || DRs->stat == bnd_uplo || DRs->stat == bnd_uld )
        && ( DRt->stat == bnd_lower || DRt->stat == bnd_uplo || DRt->stat == bnd_uld ) )
    {                          /* compare */
	lochk = ( DRs->lo_bnd == DRt->lo_bnd ) ? ck_same: ck_dif;
	loval = DRs->lo_bnd;
	if ( DRs->lo_bnd < DRt->lo_bnd )
	    loval = DRt->lo_bnd;
    }
    else
    {                          /* use the unique value */
	if ( DRs->stat == bnd_lower || DRs->stat == bnd_uplo || DRs->stat == bnd_uld )
	    loval = DRs->lo_bnd;
	else if ( DRt->stat == bnd_lower || DRt->stat == bnd_uplo || DRt->stat == bnd_uld )
	    loval = DRt->lo_bnd;
	else {
	    lo = false;        /* no lower bound */
        }
    }

/* upper bound check */

    up = true;                 /* most of the time we have an upper bound */
    upval = -99.9; /* Just for compiler warning */
    if (   ( DRs->stat == bnd_upper || DRs->stat == bnd_uplo || DRs->stat == bnd_uld )
        && ( DRt->stat == bnd_upper || DRt->stat == bnd_uplo || DRt->stat == bnd_uld ) )
    {                          /* compare */
	upchk = ( DRs->up_bnd == DRt->up_bnd ) ? ck_same: ck_dif;
	upval = DRs->up_bnd;
	if ( DRs->up_bnd > DRt->up_bnd )
	    upval = DRt->up_bnd;
    }
    else
    {                          /* use the unique value */
	if ( DRs->stat == bnd_upper || DRs->stat == bnd_uplo || DRs->stat == bnd_uld )
	    upval = DRs->up_bnd;
	else if ( DRt->stat == bnd_upper || DRt->stat == bnd_uplo || DRt->stat == bnd_uld )
	    upval = DRt->up_bnd;
	else
	    up = false;        /* no upper bound */
    }

/* compare new upper and lower bound; if incompatible then return */

    if ( upval < loval )
    {
	puts( " * Error: incompatible duplication of restraints:" );
	puts( " * Upper bound < lower bound" );
	CleanupDR_err2( count, DRs, kount, DRt );
	return( false );
    }

/* print messages */

    if ( lochk == ck_no && upchk == ck_no )
    {
	(*ndoubl)++;
        puts( "Combining bounds of:" );
	CleanupDR_lst2( count, DRs, kount, DRt );
    }
    else if (   ( lochk == ck_same || lochk == ck_no )
             && ( upchk == ck_same || upchk == ck_no ) )
    {
	(*nequal)++;
        puts( " * Repeated restraint:" );
	CleanupDR_err2( count, DRs, kount, DRt );
    }
    else
    {
	(*nequal)++;
	puts( " * Bound values differ in identical restraints (keeping the tightest bounds):" );
	CleanupDR_err2( count, DRs, kount, DRt );
    }

/* set the return values */

    *newstat = bnd_uplo;
    *newlo = loval;
    *newup = upval;

    if ( !lo )
    {
        *newlo = 0.F;
	*newstat = bnd_upper;
    }
    if ( !up )
    {
        *newup = large_BOUND;
	*newstat = bnd_lower;
    }

    return( true );
}

/****************  CleanupDR_err1  ****************/
/*
**  
*/

static void CleanupDR_err1( const int count, DRestraint *DRs )
{
    printf( " *   restraint # %i: ", count );
    WriteDR( stdout, DRs );
}

/****************  CleanupDR_lst1  ****************/
/*
**  
*/

static void CleanupDR_lst1( const int count, DRestraint *DRs )
{
    printf( "   restraint # %i: ", count );
    WriteDR( stdout, DRs );
}

/****************  CleanupDR_err2  ****************/
/*
**  
*/

static void CleanupDR_err2( const int count, DRestraint *DRs, 
                     const int kount, DRestraint *DRt )
{
    printf( " *   restraint # %i: ", count );
    WriteDR( stdout, DRs );
    printf( " *   restraint # %i: ", kount );
    WriteDR( stdout, DRt );
}

/****************  CleanupDR_lst2  ****************/
/*
**  
*/

static void CleanupDR_lst2( const int count, DRestraint *DRs, 
                     const int kount, DRestraint *DRt )
{
    printf( "   restraint # %i: ", count );
    WriteDR( stdout, DRs );
    printf( "   restraint # %i: ", kount );
    WriteDR( stdout, DRt );
}

/****************  WriteDR  ****************/
/*
**  
*/

static void WriteDR( FILE *fp, DRestraint *DRs )
{
    WriteAtom( fp, &DRs->atom_A, use_ext );
    fprintf( fp, " - ");
    WriteAtom( fp, &DRs->atom_B, use_ext );
    fprintf( fp, " : %f %f %f\n", DRs->lo_bnd, DRs->up_bnd, DRs->dist );
}

/****************  CheckBoundsValTyp  ****************/
/*
**  
*/

static Bndstat CheckBoundsValTyp( DRestraint *DRt )
{
    Bndstat  s = DRt->stat;

    if ( DRt->idx_A == -1 || DRt->idx_B == -1 )
	return( bnd_noatom );

    if ( DRt->idx_A == DRt->idx_B )
	return( bnd_sameat );

    if ( ( s == bnd_lower || s == bnd_uplo || s == bnd_uld ) && 
         DRt->lo_bnd < small_BOUND )
	return( bnd_sml_lo );

    if ( ( s == bnd_upper || s == bnd_uplo || s == bnd_uld ) && 
         DRt->up_bnd > large_BOUND )
	return( bnd_big_up );

    if ( ( s == bnd_uplo || s == bnd_uld ) && DRt->up_bnd < DRt->lo_bnd )
	return( bnd_incomp );

    if ( s == bnd_uld && !( DRt->up_bnd >= DRt->dist && DRt->lo_bnd <= DRt->dist ) )
	return( bnd_incompd );

    return( s );
}

